R Notebook

This is my script to try and analyse data from scratch - DATASET 2 -> MERSCOPE

library(Seurat)

Loading required package: SeuratObject
Loading required package: sp
The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, will retire in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.
The sp package is now running under evolution status 2
     (status 2 uses the sf package in place of rgdal)

Attaching package: ‘SeuratObject’

The following objects are masked from ‘package:base’:

    intersect, t

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(ggplot2)
library(scCustomize)

scCustomize v2.1.2
If you find the scCustomize useful please cite.
See 'samuel-marsh.github.io/scCustomize/articles/FAQ.html' for citation info.

library(readr)
library(pheatmap)
library(matrixStats)
library(spdep)

Loading required package: spData
To access larger datasets in this package, install the spDataLarge package with:
`install.packages('spDataLarge', repos='https://nowosad.github.io/drat/', type='source')`
Loading required package: sf
Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

library(geojsonR)

Read in DAta

ImageDimPlot(seurat,fov='COLON')
ImageDimPlot(seurat,fov='COLON')

# Mark the blank probes that are detected in each cell - don't think we have the other codeword categories that are present 
seurat[["Negative.Control.Codeword"]] <- CreateAssayObject(counts =data$transcripts[grepl('Blank',rownames(data$transcripts)),])

ImageFeaturePlot(seurat, "nCount_MERSCOPE", axes = T) + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

# crop only does image, not dataset
#cropped <- Crop(seurat[["COLON"]], x = c(8000, 11000), y = c(2500, 5500), coords = "plot")
#seurat[["ROIIW"]] <- cropped

# MERSCOPE CELL GETS read in as an integer but the integers are too big - need to read in cell ids as character so use tidyverse to read in set as character
global_coordinates <- data.frame(seurat@images$COLON$centroids)
global_coordinates

## This does not work - need to order correctly 
seurat$global_X <- global_coordinates$x
seurat$global_Y <- global_coordinates$y
seurat_subset <- seurat[,seurat$global_X < 11000 & seurat$global_X > 8000 & seurat$global_Y < 5500 & seurat$global_Y > 2500]

Warning: Not validating Centroids objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating Seurat objects

ImageFeaturePlot(seurat_subset, "nCount_MERSCOPE", axes = T) + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

ImageFeaturePlot(seurat_subset, "nFeature_MERSCOPE") + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

ImageFeaturePlot(seurat_subset, "volume") + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

seurat_subset[["SIZE_FILTER_LARGE"]] <- seurat_subset$volume < quantile(seurat_subset$volume, .99)
ImageDimPlot(seurat_subset, group.by="SIZE_FILTER_LARGE")

seurat_subset[["SIZE_FILTER_SMALL"]] <- seurat_subset$volume > quantile(seurat_subset$volume, .01)
ImageDimPlot(seurat_subset, group.by="SIZE_FILTER_SMALL")

```{r}

Error: attempt to use zero-length variable name

seurat_subset$TRANSCRIPT_FILTER <- seurat_subset$nCount_MERSCOPE >= 15
ImageDimPlot(seurat_subset, group.by="TRANSCRIPT_FILTER")

ImageFeaturePlot(seurat_subset, "nCount_Negative.Control.Codeword") + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

ImageFeaturePlot(seurat_subset, "nCount_Negative.Control.Probe") + scale_fill_viridis_c()

Error in `FetchData()`:
! None of the requested variables were found: nCount_Negative.Control.Probe
Backtrace:
 1. Seurat::ImageFeaturePlot(seurat_subset, "nCount_Negative.Control.Probe")
 3. SeuratObject:::FetchData.Seurat(object = object, vars = c(features, split.by[1L]), cells = cells)

ImageFeaturePlot(seurat_subset, "nFeature_Negative.Control.Codeword") + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

seurat_subset_filtered <- subset(seurat_subset, SIZE_FILTER_LARGE & SIZE_FILTER_SMALL & TRANSCRIPT_FILTER)

Warning: Not validating Centroids objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating Seurat objects

seurat_subset_filtered <- SCTransform(seurat_subset_filtered, assay = "MERSCOPE", clip.range = c(-10, 10))

Running SCTransform on assay: MERSCOPE
Running SCTransform on layer: counts
vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
Variance stabilizing transformation of count matrix of size 550 by 91161
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 549 genes, 5000 cells
Found 63 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 550 genes
Computing corrected count matrix for 550 genes
Calculating gene attributes
Wall clock passed: Time difference of 15.35141 secs
Determine variable features
Centering data matrix

  |                                                                                                     
  |                                                                                               |   0%
  |                                                                                                     
  |===============================================================================================| 100%
Getting residuals for block 1(of 19) for counts dataset
Getting residuals for block 2(of 19) for counts dataset
Getting residuals for block 3(of 19) for counts dataset
Getting residuals for block 4(of 19) for counts dataset
Getting residuals for block 5(of 19) for counts dataset
Getting residuals for block 6(of 19) for counts dataset
Getting residuals for block 7(of 19) for counts dataset
Getting residuals for block 8(of 19) for counts dataset
Getting residuals for block 9(of 19) for counts dataset
Getting residuals for block 10(of 19) for counts dataset
Getting residuals for block 11(of 19) for counts dataset
Getting residuals for block 12(of 19) for counts dataset
Getting residuals for block 13(of 19) for counts dataset
Getting residuals for block 14(of 19) for counts dataset
Getting residuals for block 15(of 19) for counts dataset
Getting residuals for block 16(of 19) for counts dataset
Getting residuals for block 17(of 19) for counts dataset
Getting residuals for block 18(of 19) for counts dataset
Getting residuals for block 19(of 19) for counts dataset
Centering data matrix

  |                                                                                                     
  |                                                                                               |   0%
  |                                                                                                     
  |===============================================================================================| 100%
Finished calculating residuals for counts
Set default assay to SCT

seurat_subset_filtered <- RunPCA(seurat_subset_filtered)

PC_ 1 
Positive:  EPCAM, LGR5, EPHB3, RGMB, SOX9, CDH1, PKM, MYC, VEGFA, MUC1 
       AXIN2, CDCA7, IDH1, MKI67, ERBB2, PROX1, ASCL2, ERBB3, SMOC2, CTNNB1 
       EPHB4, NOTCH1, PCNA, HDAC1, EPHA2, LGALS9, MET, JUN, MCM2, CEACAM1 
Negative:  COL1A1, ACTA2, FN1, COL4A1, COL5A1, PDGFRB, PECAM1, MYH11, CD248, ETS1 
       CAV1, PLVAP, ITGA1, HLA-DRA, ITGA5, DES, ENG, MMP11, ZEB1, ANGPT2 
       DUSP1, WWTR1, AKT3, ELN, COL11A1, NRP1, ADAMTS4, SFRP2, VWF, CSF1 
PC_ 2 
Positive:  SPP1, HLA-DRA, ITGB2, CYBB, FCGR3A, C1QC, CD14, CSF1R, HLA-DRB1, MRC1 
       HLA-DQA1, LYZ, PTPRC, CIITA, HLA-DPB1, ITGAX, CXCR4, FCGR2A, HLA-DMA, SERPINA1 
       ITGAM, HAVCR2, MAFB, IL4I1, TREM2, CD4, TLR2, DUSP1, CEBPB, FOS 
Negative:  FN1, COL4A1, ACTA2, COL1A1, SMOC2, PDGFRB, COL5A1, LGR5, RGMB, ITGA1 
       CD248, MMP11, PROX1, CAV1, EPHB3, ANGPT2, CDCA7, MYH11, PLVAP, NOTCH1 
       ELN, ADAMTS4, CTNNB1, INSR, ETS1, NDUFA4L2, MKI67, ASCL2, WWTR1, SOX9 
PC_ 3 
Positive:  SMOC2, MKI67, LGR5, CDCA7, PLK1, RGMB, PCNA, PROX1, BIRC5, HLA-DRA 
       CCNB1, FOXM1, MCM6, MCM2, FN1, CDK4, EPHB3, ITGB2, C1QC, CTNNB1 
       CSF1R, HLA-DRB1, MYBL2, CYBB, ASCL2, HLA-DQA1, FCGR3A, AURKB, EZH2, CD14 
Negative:  FOS, VEGFA, EPHA2, CEACAM1, JUNB, EGR1, LAMC2, SLC26A3, CDKN1A, MUC1 
       KIT, LAMB3, CXCL1, LRP1, ATF3, NFKBIA, EPCAM, LDHA, PLOD2, PPARD 
       PKM, IFNGR2, COL4A1, DUSP1, NFKB2, INSR, TCF7L2, KDR, MMP7, ITGB1 
PC_ 4 
Positive:  FN1, LGR5, SMOC2, RGMB, ERBB3, PROX1, COL1A1, MMP11, LRP1, COL5A1 
       STAT6, ERBB2, COL11A1, DES, RORC, EPHB3, SFRP2, CDKN1B, SOX9, AXIN2 
       EPCAM, TNFSF10, LRP5, ELN, VEGFA, JUN, BCL2, IGF1R, FGFR3, EPHA4 
Negative:  MKI67, PLK1, BIRC5, PCNA, FOXM1, PLVAP, CCNB1, PECAM1, PKM, KDR 
       MYBL2, VWF, MMRN2, AURKB, ENG, MCM6, AURKA, MCM2, FLT4, COL4A1 
       FLT1, ANGPT2, EZH2, PDGFB, ADAMTS4, NRP1, ETS1, BRCA1, CLEC14A, E2F1 
PC_ 5 
Positive:  PLVAP, PECAM1, INSR, VWF, MMRN2, KDR, ENG, LGR5, FLT4, SMOC2 
       FLT1, CDH5, PDGFB, ETS1, RGMB, ANGPT2, NRP1, CLEC14A, ADAMTS4, PREX2 
       COL4A1, ERBB3, PROX1, CD40, CXCR4, ACKR3, NOTCH1, ITGA5, THBD, TNFRSF4 
Negative:  MKI67, COL1A1, COL5A1, PLK1, ACTA2, FN1, PKM, BIRC5, PCNA, FOXM1 
       DES, CCNB1, EGR1, FOS, COL11A1, TGFBI, MYH11, SFRP2, TNC, ELN 
       CDKN1A, MYBL2, AURKB, AURKA, MCM6, MCM2, JUNB, CSF1, CA9, MMP11 

ElbowPlot(seurat_subset_filtered, 50)

PC_Plotting(seurat_subset_filtered, dim_number = 8)

FeaturePlot(seurat_subset_filtered, "MRC1", reduction = "pca") + scale_color_viridis_c()

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

ImageFeaturePlot(seurat_subset_filtered, "PC_1") + scale_fill_viridis_c()

Warning: No FOV associated with assay 'SCT', using global default FOVScale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

ImageFeaturePlot(seurat_subset_filtered, "MRC1", size=.5) + scale_fill_viridis_c()

Warning: No FOV associated with assay 'SCT', using global default FOVScale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

#Add annotation

ref <- readRDS("/project/shared/spatial_data_camp/datasets/SINGLE_CELL_REFERENCES/COLON_HC_5K_CELLS.RDS")
ref
DimPlot(ref)
ref <- SCTransform(ref, residual.features =rownames(seurat_subset_filtered_UMAP))
ref <- RunPCA(ref)
ref <- RunUMAP(ref, dims=1:20)
DimPlot(ref, label=T, repel=T)
ps <- AggregateExpression(ref, features = rownames(seurat_subset_filtered_UMAP), normalization.method = "LogNormalize", assays="RNA", return.seurat = T)

#Transfer labels

ImageDimPlot(seurat_subset_filtered_UMAP, size=.5)

Warning: No FOV associated with assay 'SCT', using global default FOV

FeaturePlot(seurat_subset_filtered_UMAP, "CD3E", label=T, repel=T)+ scale_color_viridis_c(direction=-1)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

top

 [1] "LGR5"     "RORC"     "RGMB"     "SMOC2"    "ERBB3"    "MKI67"    "PLK1"    
 [8] "AURKB"    "CCNB1"    "FOXM1"    "SLC26A3"  "MMP7"     "CXCL1"    "CEACAM1" 
[15] "LAMC2"    "FN1"      "CDH5"     "PROX1"    "BCL2"     "PDGFRB"   "COL5A1"  
[22] "COL1A1"   "ELN"      "DES"      "CLCA1"    "MUC2"     "SERPINA1" "Blank-46"
[29] "PTGDR2"   "SPP1"     "CSF1R"    "HLA-DRA"  "ITGB2"    "MMP12"    "PLVAP"   
[36] "VWF"      "PECAM1"   "MMRN2"    "FLT4"     "CD3E"     "TRAC"     "CD2"     
[43] "CCL5"     "EOMES"   

Clustered_DotPlot(seurat_subset_filtered_UMAP, features = top, group.by = "predicted.id", k=9)

[[1]]

[[2]]

ImageDimPlot(seurat_subset_filtered_UMAP, fov="ROI", boundaries="segmentation", border.color = "black" )

Error in ImageDimPlot(seurat_subset_filtered_UMAP, fov = "ROI", boundaries = "segmentation",  : 
  No compatible spatial coordinates present

ImageDimPlot(seurat_subset_filtered_UMAP, group.by = "predicted.id")

Warning: No FOV associated with assay 'SCT', using global default FOV

cell_colours <- c("#F8766D", "#DB8E00", "#AEA200", "#64B200", "#00BD5C", "#00C1A7", 
                  "#00BADE", "#00A6FF", "#B385FF", "#EF67EB", "#FF63B6")
names(cell_colours)  <- c("Epithelium", "Fibroblasts", "T-Cells",  "Myofibroblasts", "Macrophages", 
                          "Glia", "Endothelium", "Telocytes", "Plasma", "B-Cells", "Pericytes")
cell_colours

    Epithelium    Fibroblasts        T-Cells Myofibroblasts    Macrophages 
     "#F8766D"      "#DB8E00"      "#AEA200"      "#64B200"      "#00BD5C" 
          Glia    Endothelium      Telocytes         Plasma        B-Cells 
     "#00C1A7"      "#00BADE"      "#00A6FF"      "#B385FF"      "#EF67EB" 
     Pericytes 
     "#FF63B6" 

ImageFeaturePlot(seurat_subset_filtered_UMAP, "SMOC2") + scale_fill_viridis_c()

Warning: No FOV associated with assay 'SCT', using global default FOVScale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

markers <- FindAllMarkers(seurat_subset_filtered_UMAP, group.by="predicted.id")

explore the object

Explore QC Metrics

Filter poor quality cells

Cluster the data

Transfer the labels onto Xenium from sc RNA reference dataset - what works well and what does not

how much of it is limitations of the panel